/**
 * \file pvt_dead_oil.cpp
 * \brief
 * \author Miryanov Sergey
 * \date 07.05.2008
 */
#include "bs_pvt_stdafx.h"

#include "pvt_dead_oil.h"
#include "pvt_interpolator.h"
#include "scal_interpolate.h"

using namespace blue_sky::pvt;

namespace blue_sky
  {

  template <typename strategy_t>
  pvt_dead_oil<strategy_t>::pvt_dead_oil (bs_type_ctor_param)
  {
  }

  template <typename strategy_t>
  pvt_dead_oil<strategy_t>::pvt_dead_oil (const pvt_dead_oil &pvt)
  : bs_refcounter (pvt)
  {
    if (this != &pvt)
      {
        // TODO: LOG
        *this = pvt;
        BS_ASSERT (false && "NOT IMPL");
      }
  }


  template <typename strategy_t>
  void
  pvt_dead_oil<strategy_t>::insert_vector (const input_vector_t &vec)
  {
    const int elem_count = 3;
    BS_ASSERT (!(vec.size() % elem_count)) (vec.size ()) (elem_count);

    if (!base_t::init_dependent)
      {
        base_t::init_dependent = true;

        main_gpr_.clear ();
        main_pressure_.clear ();
        main_fvf_.clear ();
        main_visc_.clear ();
      }

    BS_ASSERT (main_gpr_.empty ());
    for (int i = 0, cnt = (int)(vec.size() / elem_count); i < cnt; ++i)
      {
        main_gpr_.push_back       (0.0);
        main_pressure_.push_back  (vec[i * elem_count + 0]);
        main_fvf_.push_back       (vec[i * elem_count + 1]);
        main_visc_.push_back      (vec[i * elem_count + 2]);
      }
  }

  template <typename strategy_t>
  int
  pvt_dead_oil<strategy_t>::build_internal (item_t atm_p, item_t min_p, item_t max_p, int n_intervals, bool is_pvto)
  {
    check_oil ();

    if (main_pressure_.empty ())
      return 0;

    index_t n_points = (index_t)main_pressure_.size ();
    if (is_pvto)
      {
        n_points = 0;
        for (index_t i = 0, cnt = (index_t)main_gpr_.size (); i < cnt; ++i)
          {
            if (main_gpr_[i] > -0.5)
              ++n_points;
          }
      }

    bool is_min         = false,      is_max    = false;
    index_t i_min       = 0,          i_max     = (index_t)main_pressure_.size () - 1;
    index_t i_next_min  = i_min + 1,  i_pre_max = i_max - 1;

    if (is_pvto)
      {
        while (main_gpr_[i_min] < 0.)
          ++i_min;

        if (min_p < main_pressure_[i_min])
          {
            is_min = true;
            n_points += (is_min != false);

            i_next_min = i_min + 1;
            while (main_gpr_[i_next_min] < 0.)
              ++i_next_min;

            if (i_next_min > (index_t)main_pressure_.size () - 1)
              {
                BS_ASSERT (i_next_min <= (index_t)main_pressure_.size () - 1) (i_next_min) (main_pressure_.size ());
                throw bs_exception ("pvt_dead_oil::build", "Error: number of Rs values in PVTO should be greater than 1");
              }
          }

        while (main_gpr_[i_max] < 0.)
          --i_max;

        if (max_p > main_pressure_[i_max])
          {
            is_max = true;
            n_points += (is_max != false);

            i_pre_max = i_max - 1;
            while (main_gpr_[i_pre_max] < 0.)
              --i_pre_max;

            if (i_pre_max < 0)
              {
                BS_ASSERT (i_pre_max >= 0) (i_pre_max);
                throw bs_exception ("pvt_dead_oil::build", "Error: number of Rs values in PVTO should be greater than 1");
              }
          }
      }
    else
      {
        is_min = min_p < main_pressure_.front ();
        is_max = max_p > main_pressure_.back ();
        n_points += (is_min != false);
        n_points += (is_max != false);
      }

    pressure_.assign (n_points, 0);
    inv_fvf_.assign (n_points, 0);
    inv_visc_.assign (n_points, 0);
    inv_visc_fvf_.assign (n_points, 0);
    gor_.assign (n_points, 0);

    if (main_pressure_.size () < 2)
      {
        if (is_pvto)
          {
            BS_ASSERT (main_pressure_.size () >= 2) (main_pressure_.size ());
            throw bs_exception ("pvt_dead_oil::build", "Error: number of rows in PVTO should be greater than 1");
          }

        for (index_t i = 0; i < n_points; ++i)
          {
            if (i == 0 && is_min)
              pressure_.front () = min_p;
            else if (i == n_points - 1 && is_max)
              pressure_.back () = max_p;
            else
              pressure_[i] = main_pressure_.front ();

            gor_[i]           = 0.0;
            inv_fvf_[i]       = 1.0 / main_fvf_.front ();
            inv_visc_[i]      = 1.0 / main_visc_.front ();
            inv_visc_fvf_[i]  = inv_fvf_[i] * inv_visc_[i];
          }
      }
    else
      {
        index_t i_rs = 0;
        for (index_t i = 0; i < n_points; ++i)
          {
            if ((i == 0 && is_min) || (i == n_points - 1 && is_max))
              {
                pressure_[i] = i == 0 ? min_p : max_p;

                item_t diff = 0;

                if (i == 0)
                  {
                    diff = (min_p - main_pressure_[i_min]) / (main_pressure_[i_next_min] - main_pressure_[i_min]);
                    inv_fvf_[i] = 1.0 / (main_fvf_[i_min] + (main_fvf_[i_next_min] - main_fvf_[i_min]) * diff);
                    inv_visc_[i] = 1.0 / (main_visc_[i_min] + (main_visc_[i_next_min] - main_visc_[i_min]) * diff);

                    if (is_pvto)
                      gor_[i] = main_gpr_[i_min] + (main_gpr_[i_next_min] - main_gpr_[i_min]) * diff;

                    if (!is_pvto || gor_[i] < 0)
                      gor_[i] = 0;
                  }
                else
                  {
                    diff = (max_p - main_pressure_[i_max]) / (main_pressure_[i_max] - main_pressure_[i_pre_max]);
                    inv_fvf_[i] = 1.0 / (main_fvf_[i_max] + (main_fvf_[i_max] - main_fvf_[i_pre_max]) * diff);
                    inv_visc_[i] = 1.0 / (main_visc_[i_max] + (main_visc_[i_max] - main_visc_[i_pre_max]) * diff);

                    if (inv_fvf_[i] < 0)
                      inv_fvf_[i] = inv_fvf_[i - 1];
                    if (inv_visc_[i] < 0)
                      inv_visc_[i] = inv_visc_[i - 1];

                    if (is_pvto)
                      gor_[i] = main_gpr_[i_max] + (main_gpr_[i_max] - main_gpr_[i_pre_max]) * diff;

                    if (!is_pvto || gor_[i] < 0)
                      gor_[i] = 0;
                  }

                inv_visc_fvf_[i] = inv_fvf_[i] * inv_visc_[i];
              }
            else
              {
                if (is_pvto)
                  {
                    while (main_gpr_[i_rs] < -0.5)
                      ++i_rs;
                  }

                pressure_[i]      = main_pressure_[i_rs];
                gor_[i]           = main_gpr_[i_rs];
                inv_fvf_[i]       = 1.0 / main_fvf_[i_rs];
                inv_visc_[i]      = 1.0 / main_visc_[i_rs];
                inv_visc_fvf_[i]  = inv_fvf_[i] * inv_visc_[i];
                ++i_rs;
              }
          }
      }

    return n_points;
  }

  template <typename strategy_t>
  void
  pvt_dead_oil<strategy_t>::build (item_t atm_p, item_t min_p, item_t max_p, int n_intervals)
  {
    build_internal (atm_p, min_p, max_p, n_intervals, false);
  }

  template <typename strategy_t>
  void
  pvt_dead_oil<strategy_t>::check_oil ()
  {
    base_t::check_common ();

    check_oil_common (main_pressure_, main_fvf_, main_visc_);
    check_gas_common (main_pressure_, main_fvf_, main_visc_);
  }

  template <typename strategy_t>
  bool
  pvt_dead_oil<strategy_t>::calc (const bool is_g, const int main_var, const item_t p, const item_t gor,
                                  item_t *inv_fvf, item_t *d_inv_fvf, item_t *gor_d_inv_fvf,
                                  item_t *inv_visc, item_t *d_inv_visc, item_t *gor_d_inv_visc,
                                  item_t *inv_visc_fvf, item_t *d_inv_visc_fvf, item_t *gor_d_inv_visc_fvf,
                                  item_t *gas_oil_ratio, item_t *d_gas_oil_ratio,
                                  const item_t drsdt /* = -1.0 */, const item_t dt /* = 0 */,
                                  const item_t old_gas_oil_ratio /* = 0  */) const
    {
      return calc_saturated_oil (is_g, main_var, p, gor, inv_fvf, d_inv_fvf, gor_d_inv_fvf,
                                 inv_visc, d_inv_visc, gor_d_inv_visc,
                                 inv_visc_fvf, d_inv_visc_fvf, gor_d_inv_visc_fvf,
                                 gas_oil_ratio, d_gas_oil_ratio, drsdt, dt, old_gas_oil_ratio);
    }

  template <typename strategy_t>
  bool
  pvt_dead_oil<strategy_t>::calc_saturated_oil (const bool is_g, const int main_var, const item_t p, const item_t gor,
      item_t *inv_fvf, item_t *d_inv_fvf, item_t *gor_d_inv_fvf,
      item_t *inv_visc, item_t *d_inv_visc, item_t *gor_d_inv_visc,
      item_t *inv_visc_fvf, item_t *d_inv_visc_fvf, item_t *gor_d_inv_visc_fvf,
      item_t *gas_oil_ratio, item_t *d_gas_oil_ratio,
      const item_t drsdt /* = -1.0 */, const item_t dt /* = 0 */,
      const item_t old_gas_oil_ratio /* = 0 */) const
    {
      size_t i = binary_search (p, pressure_, std::less <item_t> ());
      if (i == 0)
        ++i;

      if (i >= pressure_.size ())
        {
          BOSERR (section::pvt, level::error) << "pvt_dead_oil::calc_saturated_oil: index_out_of_range: " << i << " [" << p << "]" << bs_end;
          return false;
        }

      item_t dp           = pressure_[i] - pressure_[i - 1];
      item_t idp          = 1.0 / dp;
      item_t diff_p       = p - pressure_[i - 1];

      item_t d_ifvf       = (inv_fvf_[i]          - inv_fvf_[i - 1]) * idp;
      item_t ifvf         = (inv_fvf_[i - 1]      + d_ifvf * diff_p);

      item_t d_ivisc      = (inv_visc_[i]         - inv_visc_[i - 1]) * idp;
      item_t ivisc        = (inv_visc_[i - 1]     + d_ivisc * diff_p);

      item_t d_ivisc_fvf  = (inv_visc_fvf_[i]     - inv_visc_fvf_[i - 1]) * idp;
      item_t ivisc_fvf    = (inv_visc_fvf_[i - 1] + d_ivisc_fvf * diff_p);

      item_t d_gor        = (gor_[i] - gor_[i - 1]) * idp;
      item_t cur_gor      = (gor_[i - 1] + d_gor * diff_p);

      set_pvt_pointer (inv_fvf, ifvf);
      set_pvt_pointer (inv_visc, ivisc);
      set_pvt_pointer (inv_visc_fvf, ivisc_fvf);
      set_pvt_pointer (d_inv_fvf, d_ifvf);
      set_pvt_pointer (d_inv_visc, d_ivisc);
      set_pvt_pointer (d_inv_visc_fvf, d_ivisc_fvf);

      set_pvt_pointer (gas_oil_ratio, cur_gor);
      set_pvt_pointer (d_gas_oil_ratio, d_gor);

      set_pvt_pointer (gor_d_inv_fvf, 0);
      set_pvt_pointer (gor_d_inv_visc, 0);
      set_pvt_pointer (gor_d_inv_visc_fvf, 0);

      return true;
    }

  template <typename strategy_t>
  typename pvt_dead_oil<strategy_t>::item_t
  pvt_dead_oil<strategy_t>::interpolate_and_fix (item_t cell_pbub) const
    {
      size_t l = binary_search (cell_pbub, pressure_, std::less <item_t> ());
      size_t n = pressure_.size ();
      if (n == 1)
        {
          if (l == n)
            return gor_.back ();
          else if (l == 0)
            return gor_.front ();
          else
            throw bs_exception ("pvt_oil::interpolate_and_fix", "invalid index");
        }
      else if (n > 0)
        {
          if (l == n)
            l--;
          else if (l == 0)
            l++;

          item_t g2 = gor_[l - 1], g1 = gor_[l];
          item_t p2 = pressure_[l - 1], p1 = pressure_[l];
          return g2 + (cell_pbub - p2) * (g1 - g2)	/ (p1 - p2);
        }
      else
        {
          bs_throw_exception (boost::format ("Unexpected error: n (%lu) out of range") % n);
        }
    }

  template <typename strategy_t>
  typename pvt_dead_oil<strategy_t>::item_t
  pvt_dead_oil<strategy_t>::get_gor_for_pressure (item_t pressure_data) const
    {
      if (pressure_data < pressure_.front ())
        {
          throw bs_exception ("pvt_oil::get_gor_for_pressure", "Invalid pressure data value");
        }

      index_t i_rs = (index_t) (((pressure_data) - pressure_.front ()) / this->get_p_step ());

      return gor_[i_rs] + (gor_[i_rs + 1] - gor_[i_rs]) / this->get_p_step () * (pressure_data - pressure_[i_rs]);
    }

  template <typename strategy_t>
  void
  pvt_dead_oil <strategy_t>::print () const
  {
    BS_ASSERT (pressure_.size () == inv_fvf_.size ());
    BS_ASSERT (inv_fvf_.size ()  == inv_visc_.size ());
    BS_ASSERT (inv_visc_.size () == inv_visc_fvf_.size ());
    
    BOSOUT (section::pvt, level::medium) << bs_end;
    BOSOUT (section::pvt, level::medium) << "*************************************** PVDO ****************************************" << bs_end;
    BOSOUT (section::pvt, level::medium) 
      << boost::format ("*%9s%13s%12s%12s%12s%12s%12s%12s%2s\n")
      % "P"
      % "RS"
      % "FVF"
      % "1/FVF"
      % "VISC"
      % "1/VISC"
      % "VISC*FVF"
      % "1/FVF/VISC"
      % "*"
      << bs_end;
    BOSOUT (section::pvt, level::medium) << "*************************************************************************************" << bs_end;
    for (size_t j = 0, jcnt = pressure_.size ();j < jcnt; ++j)
      {
        BOSOUT (section::pvt, level::medium)
          << boost::format (
            "*%9.6f\t\t%9.6f\t\t%9.6f\t\t%9.6f\t\t%9.6f\t\t%9.6f\t\t%9.6f\t\t%9.6f%2s"
          ) 
          % (pressure_[j])
          % (gor_[j])
          % (1.0 / inv_fvf_[j])
          % (inv_fvf_[j])
          % (1.0 / inv_visc_[j])
          % (inv_visc_[j])
          % (1.0 / (inv_visc_fvf_[j]))
          % (inv_visc_fvf_[j])
          % "*"
          << bs_end;
      }
    BOSOUT (section::pvt, level::medium) << "*************************************************************************************" << bs_end;
  }

  template class pvt_dead_oil <base_strategy_fi>;
  template class pvt_dead_oil <base_strategy_di>;
  template class pvt_dead_oil <base_strategy_mixi>;
} // namespace blue_sky

